#installation de package nécéssaire en local
#install.packages("BiocManager")
#BiocManager::install("phyloseq")
#install.packages("phyloseq")
#install.packages("remotes")
#remotes::install_github("mahendra-mariadassou/phyloseq-extended", ref = "dev")
#install.packages("gridExtra")
#Installation des librairies
library(tidyverse)
library(phyloseq)
library(phyloseq.extended)
library(ggplot2) ## for graphs necessary ???
library(ape) ## for tree manipulation
library(vegan) ## for community ecology analyses
library(scales)
library(grid)
library(gridExtra) ## for design graph
library(knitr)
library(kableExtra)
#chemin <- "~/save/Milint/Resultat/"
load("BIOM.Rdata")
#load("Data_Biom.Rdata")Le BIOM a été créé dans le script01. Il est cependant “brut”.
Somme de tous les reads et on applique un filtre de 10-5 afin de s’assurer que ce ne sont pas des artefacts Dans notre cas permet de filtrer les taxa avec des valeurs supérieures à 180 reads par OTU
filtre <- sum(taxa_sums(Biom))*1e-5
filter_biom <- prune_taxa(taxa_sums(Biom)>filtre,Biom)
save(Biom, filter_biom, file = "Data_Biom.Rdata")
[1] 5042 240312
[1] 5003 239715
# moyenne
moy <- sum(taxa_sums(Biom))/846
moyf <- sum(taxa_sums(filter_biom))/846
mi <- min(taxa_sums(Biom))
ma <- max(taxa_sums(Biom))
mif <- min(taxa_sums(filter_biom))
maf <- max(taxa_sums(filter_biom))
med <- median(taxa_sums(Biom))
medf <- median(taxa_sums(filter_biom))
df <- data.frame(
échantillon = c(846,846),
mean = c(moy,moyf),
median = c(med,medf),
min = c(mi,mif),
max = c(ma, maf)
)
rownames(df) <- c("Biom","Filter_Biom")
kable(df)| échantillon | mean | median | min | max | |
|---|---|---|---|---|---|
| Biom | 846 | 21339.08 | 5 | 0 | 2709474 |
| Filter_Biom | 846 | 21207.16 | 1092 | 181 | 2709474 |
Le niveau d’abondance semble hétérogène. Il faudra donc prévoir un procédé de rarefaction
plot_bar(Biom)+
ggtitle("Abondance par échantillon sur données totales")+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())plot_bar(filter_biom)+
ggtitle("Abondance par échantillon sur données filtrées")+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())Dorénavant nous ne travaillerons que sur données filtrées Ne sert à rien d’aller au niveau Phylum ou infra car trop d’échantillons donc aucune couleur n’est visible Nous remarquons que le niveau d’abondance par échantillon est hétérogène de 25 000 en moyenne à 250 000.
Analyses effectuées mais retirées : J’avais fait des niveaux d’abondance par échantillon en fonction de différentes variables (sexe, satut tabagique…) mais ne sert à rien car ça équivaut à voir la profondeur et ce n’est pas biologiquement informatif donc enlevé.
| Var1 | Freq |
|---|---|
| Actinobacteria | 22 |
| Bacteroidetes | 223 |
| Firmicutes | 686 |
| Fusobacteria | 2 |
| Proteobacteria | 45 |
| Synergistetes | 9 |
| Tenericutes | 5 |
| Verrucomicrobia | 5 |
Les Firmicutes sont les bactéries majoritaires. Elles sont les plus abondantes suivies des Bacteroidetes
#Abondance par phylum
phy <- tax_glom(filter_biom, taxrank = "Phylum")
plotdata<-psmelt(phy)
ggplot(plotdata,aes(y=Abundance))+
geom_boxplot(aes(group="Phylum", fill="Phylum"))+
facet_wrap(~Phylum, scales="free_y",nrow=2)plot_composition(filter_biom, "Phylum", "Firmicutes", "Family", fill = "Family")+
facet_grid(scales = "free_x") +
theme(axis.text.x = element_blank())Pas de différence flagrante. Les échantillons semblent assez hommogène dans la composition par famille. Sur les données au niveau du genre rien ne se dégage ainsi qu’en fonction du sexe, age, tabac… Les échantillons semblent peu différents entre eux. Nous risquons donc d’avoir une faible voir une absence de diversité béta
Au vu des différents graphs, il y a une grande diversité de bactéries et nous ne voyons pas de genre majoritaire quelque soit la variable. Pour l’ordre Clostrdiales est majoritaires quelque soit la variable prise (exemple ci dessous)
Composition en Firmicutes au niveau de l’ordre en fonction de l’âge
plot_composition(filter_biom, "Phylum", "Firmicutes","Order", fill = "Order") +
facet_grid(~AGE, scales = "free_x", space = "free_x")+
theme(axis.text.x = element_blank())Je pense que les graphs suivants ne sont d’aucune utilité . J’ai fait les mêmes sur les femmes mais trop d’échantillons
Diversité intrinsèque à chaque échantillon
Exploration de l’impact de chaque covariable sur la diversité-alpha
Shannon, Cet indice représente à la fois le nombre d’espèces d’un milieu mais aussi la répartition des effectifs individuels au sein des espèces présentes. Simpson mesure de régularité cad mesure la probabilité que deux individus pris au hasard appartiennent à la même espèce Chao1 tient davantage compte des espèces peu abondantes nécéssite de conserver singleton donc pas pris
plot_richness(filter_rare, x = "AGE", measures = c("Observed", "Shannon", "InvSimpson")) +
geom_boxplot(aes(group = AGE)) ## add one boxplot per AGE valuesAge diffère significativement en terme d’observation d’OTU
div_data <- cbind(estimate_richness(filter_biom), ## diversity indices
sample_data(filter_biom) ## covariates
)Mais également en terme de diversité en nombre d’espèce
plot_richness(filter_rare, x = "SEX", measures = c("Observed", "Shannon", "InvSimpson")) +
geom_boxplot(aes(group = SEX)) ## add one boxplot per SEX valuesPas de différence signicative pour le nbre d’OTU mais différence significative pour le nbre d’espèces
plot_richness(filter_rare, x = "tabac", measures = c("Observed", "Shannon", "InvSimpson")) +
geom_boxplot(aes(group = tabac)) ## add one boxplot per AGE valuesObservation significative ainsi que Shannon. InvSimpson 0.08
div_data <- cbind(estimate_richness(filter_rare), ## diversity indices
sample_data(filter_rare)) ## covariates
model <- aov(Observed ~ tabac+age+SEX+BMI, data = div_data)
anova(model), ,) SEX nscorrespond à la différence de diversité des espèces entre plusieurs milieux Diversité entre échantillons : indice de dissimilarité (Bray et Curtis , de Jacard) Phylogénie (Unifrac) ### Nous n’avons pas d’arbre phylogénétique donc pas de distance unifrac et wunifrac
Conclusion : Aucune tendance se dégage. Il n’y a pas de différence de diversité d’espèces entre les différent échantillons pour les catégories (age, sexe, tabac, bmi). Impossible de voir pour l’activité physique. La plage semble trop importante pour une coloration. Voir pour mettre des catégories
dist.jac <- distance(filter_rare, method = "jaccard")
dist.bc <- distance(filter_rare, method = "bray")
#dist.uf <- distance(xxx, method = "unifrac")
#dist.wuf <- distance(xxx, method = "wunifrac")Représentation MDS (pour « Metric MultiDimensional Scaling » soit analyse multidimensionnelle métrique) La NMDS ne converge pas et NMDS déforme l’espace pour faire apparaitre des groupes éventuels. Du coup les distances apparentes ne sont pasfidèles aux distances réelles.
my_palette <- c('red','yellow','purple','green','orange')
p.jac <- plot_ordination(filter_rare,
ordinate(filter_rare, method = "MDS", distance = dist.jac),
color = "age") +
ggtitle("Jaccard")+ scale_color_manual(values = my_palette)
p.bc <- plot_ordination(filter_rare,
ordinate(filter_rare, method = "MDS", distance = dist.bc),
color = "age") +
ggtitle("Bray-Curtis")+ scale_color_manual(values = my_palette)
gridExtra::grid.arrange(p.jac, p.bc, ncol = 2)Quelque soit la variable utilisée, aucun structuration n’est identifiée. Les 2 premiers axes ne capturent que que 22% de la diversité avec la distance de Bray-Curtis.
/ /
Risque de ne rien voir car déjà rien ne ressort sans contrainte Changement de l’argument method en CAP (Constrained Analysis of Proximities) et indiqué quelle covariable peut expliquer cette diversité
ord <- ordinate(filter_rare, method = "CAP", distance = dist.bc, formula = ~ age)
p12 <- plot_ordination(filter_rare, ord, color = "age", axes = c(1:2))
p23 <- plot_ordination(filter_rare, ord, color = "age", axes = c(2:3))
gridExtra::grid.arrange(p12, p23, nrow = 1)/
Aucune clusterisation des échantillons n’est visible Si présence d’arbre phylo possibilité de faire avec en utilisant dist=“unifrac”
par(mar = c(1, 0, 2, 0))
plot_clust(filter_rare, dist = "bray", method = "ward.D2", color = "SEX",
title = "Clustering of samples (Bray-Curtis + Ward)\nsamples colored by sexe")/ /
metadata <- sample_data(filter_rare) %>% as("data.frame")
model <- vegan::adonis(dist.bc ~ age + SEX + tabac + BMI + APhysGlobHParSem , data = metadata, permutations = 999)
model
Call:
vegan::adonis(formula = dist.bc ~ age + SEX + tabac + BMI + APhysGlobHParSem, data = metadata, permutations = 999)
Permutation: free
Number of permutations: 999
Terms added sequentially (first to last)
Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
age 4 1.732 0.43289 2.0403 0.00954 0.001 ***
SEX 1 1.273 1.27315 6.0007 0.00701 0.001 ***
tabac 2 0.672 0.33604 1.5839 0.00370 0.013 *
BMI 1 0.318 0.31767 1.4973 0.00175 0.065 .
APhysGlobHParSem 1 0.198 0.19769 0.9318 0.00109 0.552
Residuals 836 177.370 0.21216 0.97691
Total 845 181.562 1.00000
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Toutes les variables, à l’exception de l’activité physique, sont significatives. Cependant elles n’expliquent quasiment pas de variabilités < 1%
/ / ## Heatmap
Ne voyant pas de cluster, il n’est pas nécéssaire de faire une heatmap d’une variable en considérant les distances de Bray-Curtis.
plot_heatmap(filter_rare, low = "yellow", high = "red", na.value = "white") +
facet_grid(~age, scales = "free_x", space = "free_x")+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank()) Quelque soit la tranche d’âge, il y a beaucoup de diversité. Ce sont les mêmes OTU qui sont présents dans tous les échantillons quelque soit la catégorie d’âge.
# Possibilité de zoomer sur une zone les 50 taxa les plus élevés
top_50_taxa <- taxa_sums(filter_rare) %>% sort(decreasing = TRUE) %>% names() %>% head(n = 50)
p <- plot_heatmap(prune_taxa(top_50_taxa, filter_rare),
low = "yellow", high = "red", na.value = "white") +
facet_grid(~age, scales = "free_x", space = "free_x")+
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())
plot(p)
Objectif : Voir si l’abondance de certaines espèces diffèrent entre les groupes
knitr::opts_chunk$set(fig.height = 5)
cds <- phyloseq_to_deseq2(filter_biom, ~SEX)
dds <- DESeq2::DESeq(cds, sfType = "poscounts")
results <- DESeq2::results(dds) %>% as.data.frame()
##### On ne peut prendre que les OTU dont le fold change >0.05% on peut être plus drastique
da_otus <- results %>% as_tibble(rownames = "OTU") %>%
filter(padj < 0.05) %>%
arrange(log2FoldChange) %>%
pull(OTU)
length(da_otus)[1] 60
##### Récupères ces OTU significatifs et on les ordonne dans une heatmap
p <- plot_heatmap(prune_taxa(da_otus, filter_biom), ## keep only da otus...
taxa.order = da_otus ## ordered according to fold-change
) +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())+
facet_grid(~SEX, scales = "free_x", space = "free_x") +
scale_fill_gradient2(low = "#1a9850", mid = "#ffffbf", high = "#d73027",
na.value = "white", trans = log_trans(4),
midpoint = log(100, base = 4))
plot(p)cds <- phyloseq_to_deseq2(filter_biom, ~age)
dds <- DESeq2::DESeq(cds, sfType = "poscounts")
results <- DESeq2::results(dds) %>% as.data.frame()
##### On ne peut prendre que les OTU dont le fold change >0.05% on peut être plus drastique
da_otus <- results %>% as_tibble(rownames = "OTU") %>%
filter(padj < 0.05) %>%
arrange(log2FoldChange) %>%
pull(OTU)
length(da_otus)[1] 57
##### Récupères ces OTU significatifs et on les ordonne dans une heatmap
p <- plot_heatmap(prune_taxa(da_otus, filter_biom), ## keep only da otus...
taxa.order = da_otus ## ordered according to fold-change
) +
theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank())+
facet_grid(~age, scales = "free_x", space = "free_x") +
scale_fill_gradient2(low = "#1a9850", mid = "#ffffbf", high = "#d73027",
na.value = "white", trans = log_trans(4),
midpoint = log(100, base = 4))
plotly::ggplotly(p)